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A major goal of neuroscience, statistical physics and nonlinear dynamics is to understand how brain function 
arises from the collective dynamics of networks of spiking neurons. This challenge has been chiefly addressed 
through large-scale numerical simulations. Alternatively, researchers have formulated mean-field theories to 
gain insight into macroscopic states of large neuronal networks in terms of the collective firing activity of the 
neurons, or the firing rate. However, these theories have not succeeded in establishing an exact correspondence 
between the firing rate of the network and the underlying microscopic state of the spiking neurons. This has 
largely constrained the range of applicability of such macroscopic descriptions, particularly when trying to de¬ 
scribe neuronal synchronization. Here we provide the derivation of a set of exact macroscopic equations for 
a network of spiking neurons. Our results reveal that the spike generation mechanism of individual neurons 
introduces an effective coupling between two biophysically relevant macroscopic quantities, the firing rate and 
the mean membrane potential, which together govern the evolution of the neuronal network. The resulting equa¬ 
tions exactly describe all possible macroscopic dynamical states of the network, including states of synchronous 
spiking activity. Finally we show that the firing rate description is related, via a conformal map, with a low¬ 
dimensional description in terms of the Kuramoto order parameter, called Ott-Antonsen theory. We anticipate 
our results will be an important tool in investigating how large networks of spiking neurons self-organize in time 
to process and encode information in the brain. 

PACS numbers: 87.19.1j 05.45.Xt 87.10.-Ed 05.65.+b 


Processing and coding of information in the brain necessar¬ 
ily imply the coordinated activity of large ensembles of neu¬ 
rons. Within sensory regions of the cortex, many cells show 
similar responses to a given stimulus, indicating a high de¬ 
gree of neuronal redundancy at the local level. This suggests 
that information is encoded in the population response and 
hence can be captured via macroscopic measures of the net¬ 
work activity [1]. Moreover, the collective behavior of large 
neuronal networks is particularly relevant given that current 
brain measurement techniques, such as electroencephalogra¬ 
phy (EEG) or functional magnetic resonance imaging (fMRI), 
provide data which is necessarily averaged over the activity of 
a large number of neurons. 

The macroscopic dynamics of neuronal ensembles has been 
extensively studied through computational models of large 
networks of recurrently coupled spiking neurons, including 
Hodgkin-Huxley-type conductance-based neurons [2] as well 
as simplified neuron models, see e.g. [3-5], In parallel, re¬ 
searchers have sought to develop statistical descriptions of 
neuronal networks, mainly in terms of a macroscopic observ¬ 
able that measures the mean rate at which neurons emit spikes, 
the firing rate [6-22]. These descriptions, called firing-rate 
equations (FREs), have been proven to be extremely useful 
in understanding general computational principles underlying 
functions such as memory [23, 24], visual processing [25-27], 
motor control [28] or decision making [29]. 

Despite these efforts, to date there is no exact theory link¬ 
ing the dynamics of a large network of spiking neurons with 
that of the firing rate. Specifically, current macroscopic de¬ 
scriptions do not offer a precise correspondence between the 
microscopic dynamics of the individual neurons, e.g. individ¬ 
ual membrane potentials, and the firing rate dynamics of the 


neuronal network. 

Indeed, FREs are generally derived through heuristic ar¬ 
guments which rely on the underlying spiking activity of the 
neurons being asynchronous and hence uncorrelated. As such, 
firing rate descriptions are not sufficient to describe network 
states involving some degree of spike synchronization. Syn¬ 
chronization is, however, an ubiquitous phenomenon in the 
brain, and its potential role in neuronal computation is the 
subject of intense research [14, 30-36]. Hence, the lack of 
firing rate descriptions for synchronous states limits the range 
of applicability of mean-field theories to investigate neuronal 
dynamics. 

Here we propose a method to derive the FREs for networks 
of heterogeneous, all-to-all coupled quadratic integrate-and- 
fire (QIF) neurons, which is exact in the thermodynamic limit, 
i.e. for large numbers of neurons. We consider an ansatz for 
the distribution of the neurons’ membrane potentials that we 
denominate the Lorentzian ansatz (LA). The LA solves the 
corresponding continuity equation exactly, making the system 
amenable to theoretical analysis. Specifically, for particular 
distributions of the heterogeneity, the LA yields a system of 
two ordinary differential equations for the firing rate and mean 
membrane potential of the neuronal population. These equa¬ 
tions fully describe the macroscopic states of the network — 
including synchronized states—, and represent the first exam¬ 
ple of an exact firing-rate description for a network of recur¬ 
rently connected spiking neurons. We finally show how the 
LA transforms, via a conformal mapping, into the so-called 
Ott-Antonsen ansatz (OA) that is extensively used to inves¬ 
tigate the low-dimensional dynamics of large populations of 
phase oscillators in terms of the Kuramoto order parameter 
[37], 
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I. MODEL DESCRIPTION 

Hodgkin-Huxley-type neuronal models can be broadly 
classified into two classes, according to the nature of their 
transition to spiking in response to an injected current [38, 39]. 
Neuronal models with so-called Class I excitability, generate 
action potentials with arbitrarily low frequency, depending on 
the strength of the applied current. This occurs when a resting 
state disappears through a saddle-node bifurcation. In con¬ 
trast, in neurons with Class II excitability the action potentials 
are generated with a finite frequency. This occurs when the 
resting state loses stability via a Hopf bifurcation. The QIF 
neuron is the canonical model for Class I neurons, and thus 
generically describes their dynamics near the spiking thresh¬ 
old [5, 40,41], Our aim here is to derive the FREs correspond¬ 
ing to a heterogeneous all-to-all coupled population of N QIF 
neurons. The correspondence is exact in the thermodynamic 
limit, i.e. when N —> oo (this convergence has been recently 
studied in [42]). 

The microscopic state of the population of QIF neurons is 
given by the membrane potentials which obey 

the following ordinary differential equations [5]: 

Vj = Vf + Ij , if Vj > V p , then Vj^V r . (1) 

Here the overdot denotes the time derivative, and Ij represents 
an input current. Each time a neuron’s membrane potential Vj 
reaches the peak value V p , the neuron emits a spike and its 
voltage is reset to the value V r . In our analysis we consider 
the limit V p = —V r — > oo. This resetting rule captures the 
spike reset as well as the refractoriness of the neuron. Without 
loss of generality we have rescaled the time and the voltage in 
( 1 ) to absorb any coefficients which would have appeared in 
the first two terms. The form for the input currents is: 


Ij — Vj + Js(t) + I(t), (2) 

where the external input has a heterogeneous, quenched com¬ 
ponent ijj as well as a common time-varying component /(f), 
and the recurrent input is the synaptic weight J times the mean 
synaptic activation s(f), which is written: 


s(t) 


1 N 

J = 1 k\tj <t 



dt'a T (t — 


tj)- 


(3) 


Here, tj is the time of the fcth spike of jth neuron, S(t) is the 
Dirac delta function, and a T (t) is the normalized synaptic ac¬ 
tivation caused by a single pre-synaptic spike with time scale 
r, e.g. a T (t) = e~ t/T /r. 


A. Continuous formulation 


a continuous random variable distributed according to a prob¬ 
ability distribution function //(//]. The total voltage density at 
time t is then given by f^° p(V\rj, t) g{rj) dp. 

The conservation of the number of neurons leads to the fol¬ 
lowing continuity equation: 

dtp + dv [(' V 2 + v + Js + I)p\ = 0, (4) 

where we have explicitly included the velocity given by ( 1 ) 
and ( 2 ). 


II. RESULTS 


The continuity equation (4) without temporal forcing 
I(t) = 0 has a trivial stationary solution. For each value 
of 77 , this solution has the form of a Lorentzian function: 
Po(V \v) oc ( V 2 + v+ Js ) _1 . Physically, the Lorentzian den¬ 
sity means that firing neurons with the same ?/ value will be 
scattered with a density inversely proportional to their speed 
( 1 ), i.e. they will accumulate at slow places and thin out at 
fast places on the V axis. In addition, for those 77 values cor¬ 
responding to quiescent neurons, the density po collapses at 
the rest state in the form of a Dirac delta function. 

Next we assume that, independently of the initial condition, 
solutions of (4) generically converge to a Lorentzian-shaped 
function, so that all relevant dynamics occur inside that lower 
dimensional space. This fact is mathematically expressed by 
the following Lorentzian Ansatz (LA) for the conditional den¬ 
sity functions: 


p{v\v,t) 


I _ x(y,t) _ 

71 [v - y{v,t)] 2 + x(i ?,£) 2 ’ 


(5) 


which is a Lorentzian function with time-dependent half¬ 
width x(v,t) and center at 7 /( 77 , t). In the following we as¬ 
sume that the LA (5) completely describes the macroscopic 
dynamics of the network of QIF neurons and postpone the 
mathematical justification of its validity to Sec. HE. 


A. Macroscopic observables: Firing rate and mean membrane 
potential 

The half-width x(v, t) of the LA has a particularly simple 
relation with the firing rate of the neuronal population (i.e. the 
number of spikes per unit time). Indeed, the firing rate for 
each 77 value at time t, r(ij, t ), can be computed by noting that 
neurons fire at a rate given by the probability flux at infinity: 
r(?7, t) = p(V -A oo|?7, t)V ( V —j oo|t/, t). The limit V —> 00 
on the right hand side of this equation can be evaluated within 
the LA, and gives the simple identity 

x(v,t) = nr(v,t). (6) 


In the thermodynamic limit N —> 00 , we drop the indices 
in Eqs. (1), (2) and denote p(V\v, t)dV as the fraction of neu¬ 
rons with membrane potentials between V and V + dV, and 
parameter ?/ at time t. Accordingly, parameter 77 becomes now 


The (total) firing rate r(t) is then 

r(t) = -f x(v,t)g(v)dv- (7) 

^ J — OO 
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FIG. 1. (color online). Analysis of the steady states of FREs (12). (a) Phase diagram: In the wedge-shaped cyan-shaded region there is 
bistability between a high and a low activity state. The boundary of the bistability region is the locus of a saddle-node bifurcation which is 
exactly obtained in parametric form: (fj, J)sn = [—n 2 r 2 — 3A 2 /(27rr) 2 , 2n 2 r + A 2 /(27r 2 r 3 )]. To the right of the dashed line, defined by 
fjf = —[J/(27r)] 2 — (VA/J) 2 , there is a stable focus (shaded regions), (b) r-fj and (c) v-y bifurcation diagrams for J/ A 1 ' 2 = 15. Square 
symbols: Results obtained from numerical simulations of QIF neurons (see Appendix A for details), (d) Phase portrait of the system in the 
bistable region ( 77 /A = —5, J / A 1 / 2 = 15, triangle in panel (a)) with three fixed points: a stable focus (with its basin of attraction shaded), a 
stable node, and a saddle point. 


Additionally, the quantity 7 /( 77 , t) is, for each value of 77 , the 
mean of the membrane potential: 

/ OO 

p(V\ri,t)V dV. ( 8 ) 

-OO 

Here we take the Cauchy principal value, defined as 
P.V. f(x)dx = liniR-j .00 J R r f(x)dx, to avoid the oth¬ 
erwise ill-defined integral. The mean membrane potential is 
then 


Note that this distribution accounts for the quenched variabil¬ 
ity in the external inputs. The fact that it is Lorentzian is 
unrelated to the LA for the density of membrane potentials. 
Assuming (11) the integrals in (7) and (9) can be evaluated 
closing the integral contour in the complex //-plane and using 
the residue theorem[43]. Notably, the firing rate and the mean 
membrane potential depend only on the value of w at the pole 
of g(rj) in the lower half 77 -plane: 

7 rr(f) + iv(t) = w{fj — zA, t). 


/ OO 

y(v,t)g(v)drj. 

-00 


(9) 


B. Firing-rate equations 


Substituting the LA (5) into the continuity equation (4), we 
find that, for each value of 77 , variables x and y must obey two 
coupled equations which can be written in complex form as 

d t w(rj, t) = i [77 + Js(t) - w{rj, t) 2 + /(f)] . (10) 


where 77 /( 77 > t) — X {V: t) + f)- Closing this equation re¬ 

quires expressing the mean synaptic activation s(t ) as a func¬ 
tion of 77 /( 77 , t). The simplest choice is to take the limit of in¬ 
finitely fast synapses, r —> 0 in (3), so that we get an equality 
with the firing rate: s{t) = r(f). This allows for the system of 
QIF neurons (l)-(3) to be exactly described by Eqs. (10) and 
(7); an infinite set of integro-differential equations if < 7 ( 77 ) is a 
continuous distribution. 

Equation (10) is useful for general distributions 77 ( 77 ) (see 
Appendix B), but a particularly sharp reduction in dimension¬ 
ality is achieved if 77 is distributed according to a Lorentzian 
distribution of half-width A centered at fj: 


g(v) 


1 A 

7T (77 — fj) 2 + A 2 


(ID 


As a result, we only need to evaluate (10) at 77 = 77 — iA, and 
thereby obtain a system of FREs composed of two ordinary 
differential equations 

f = A/ir + 2rv , (12a) 

v = v 2 + fj + Jr + I{t) — 7 r 2 r 2 . (12b) 

This nonlinear system describes the macroscopic dynamics of 
the population of QIF neurons in terms of the population firing 
rate r and mean membrane potential v. 

It is enlightening to compare mean-field model (12) with 
the equations of the spiking neurons. Note that Eq. (12b) re¬ 
sembles equations (1) and (2) for the individual QIF neuron, 
but without spike resetting. Indeed, the macroscopic firing- 
rate variable r enters as a negative feedback term in Eq. (12b), 
and impedes the explosive growth of the mean membrane po¬ 
tential v. 

This negative feedback, combined with the coupling term 
on the right hand side of (12a), describe the effective interac¬ 
tion between the firing rate and mean membrane potential at 
the network level. Therefore, the FREs (12) describe the ef¬ 
fect of the single-cell spike generation and reset mechanism at 
the network level. 

In the following we examine the dynamics described by 
Eq. (12), and show that they fully reproduce the macro¬ 
scopic dynamics of the network of QIF neurons, even during 
episodes of strong spike synchrony. 
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FIG. 2. (color online). The transient dynamics of an ensemble of QIF model neurons (l)-(2) are exactly described by the FREs (12). The 
instantaneous firing rate (Panels a,b) and mean membrane potential (Panels c,d) of the QIF neurons and the FREs are depicted in black and 
orange, respectively. Panels (e,f) show the raster plots of 300 randomly selected QIF neurons of the N = 10 4 in the ensemble. Left Panels 
(a,c,e,g): At time t = 0, a current I = 3 is applied to all neurons, and set to zero again at t = 30; stimulus I(t) shown in panel (g). Right 
panels (b,d,f,h): At time t = 0 a sinusoidal current is applied to all neurons /(f) = Io sin(u;f), with Jo = 3, uj = 7t/20; stimulus /(f) shown 
in panel (h). Parameters: J = 15, fj = —5, A = 1. 


C. Analysis of the firing-rate equations 

To begin with the analysis of Eq. (12), we first note that in 
the absence of forcing, /(f) = 0, the only attractors of (12) are 
fixed points. Figure 1(a) shows a phase diagram of the system 
as a function of the mean external drive fj and synaptic weight 
J, both normalized by the width of the input distribution [44]. 
There are three qualitatively distinct regions of the phase dia¬ 
gram: 1 - A single stable node corresponding to a low-activity 
state (white), 2 - A single stable focus (spiral) generally cor¬ 
responding to a high-activity state (gray), and 3 - A region of 
bistability between low and high firing rate (cyan; see a phase 
portrait of this region in Fig. 1(d)). Comparison of a sample 
bifurcation diagram of the fixed points from numerical simu¬ 
lation of networks of QIF neurons with that obtained from the 
FREs (12) shows an excellent correspondence, see Fig. l(b,c). 

A similar phase diagram can be readily reproduced by tra¬ 
ditional heuristic firing-rate models, with one significant qual¬ 
itative difference: the presence of a stable focus —and hence 
damped oscillations. Specifically, in the gray region of the 
phase diagram in Fig. 1(a), the system undergoes oscillatory 
decay to the stable fixed point. This oscillatory decay occurs 
as well for the high-activity state over a large extent of the 
region of bistability (cyan), see e.g. Fig. 1(d). 

The presence of damped oscillations at the macroscopic 
level reflects the transitory synchronous firing of a fraction 
of the neurons in the ensemble. This behavior is common in 
spiking neuron models with weak noise, and is not captured 
by traditional firing rate models (see e.g. [22]). 


D. Analysis of the firing-rate equations: Non-stationary inputs 

To show that the FREs (12) fully describe the macroscopic 
response of the population of QIF neurons to time-varying 


stimuli (up to finite-size effects), we consider two types of 
stimulus I(t)\ 1 - a step function and 2 - a sinusoidal forcing. 
In both cases we simulate the full system of QIF neurons and 
the FREs (12). 

Figure 2 shows the system’s response to the two different 
inputs. In both cases the system is initially (t < 0) in a 
bistable regime and set in the low activity state, with parame¬ 
ters corresponding to those of Fig. 1(d). Left panels of Fig. 2 
show the response of the system to a step current, applied at 
t = 0. The applied current is such that the system abandons 
the bistable region —see Fig. 1(a)— and approaches the high 
activity state, which is a stable focus. This is clearly reflected 
in the time series r(t), v(t), where the rate equations (12) ex¬ 
actly predict the damped oscillations exhibited by the mean- 
field of the QIF neurons. The raster plot in panel (e) shows 
the presence of the oscillations, which is due to the transitory 
synchronous firing of a large fraction of neurons in the popula¬ 
tion. Finally, at t = 30, the current is removed and the system 
converges —again, showing damped oscillations— to the new 
location of the (focus) fixed point, which clearly coexists with 
the stable node where it was originally placed (t < 0). 

The right hand panels of Fig. 2 show the response of the 
model to a periodic current, which drives the system from one 
side of the bistable region to the other. As a result, we observe 
periodic bursting behavior when the system visits the stable 
focus region of the phase diagram Fig. 1(a). 

To further illustrate the potential of the FREs (12) to pre¬ 
dict and investigate complex dynamics in ensembles of spik¬ 
ing neurons, we present here a simple situation where the sys¬ 
tem of QIF neurons exhibits macroscopic chaos. This is ob¬ 
served by increasing the frequency uj of the sinusoidal driving, 
so that the system cannot trivially follow the stable fixed point 
at each cycle of the applied current. 

Figure 3(a) shows a phase diagram obtained using Eq. (12). 
The shaded regions indicate parameter values where the rate 
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FIG. 3. (color online). Firing rate model predicts the existence of 
macroscopic chaos in a ensemble QIF neurons (1) with an injected 
periodic current /(f) = Iq sin(aif). (a) Phase diagram with the re¬ 
gions where chaos is found, obtained using the rate model (12). In the 
black-shaded region, there is only one chaotic attractor, whereas in 
the cyan region the chaotic attractor coexists with a periodic attractor. 
Dotted lines correspond to the bistability boundary of Fig. 1(a), de¬ 
picted to facilitate comparison, (b) Chaotic trajectory obtained sim¬ 
ulating the FREs (12); the Lyapunov exponent is A = 0.183 .... (c) 
Chaotic trajectory obtained from the QIF neurons —parameters cor¬ 
responding to the (yellow) triangle symbol in panel (a), (d) Time 
series for the rate model (orange) and QIF model (black), (e) Raster 
plot corresponding to 300 randomly selected neurons, (f) Same raster 
plot as in (e), ordering neurons according to their intrinsic current r/k ■ 
Parameters: Jo = 3, oj = 7r; and ij = —2.5, J = 10.5 (triangle sym¬ 
bol in (a)). 


model has either a single chaotic attractor (in black) or a 
chaotic attractor coexisting with a periodic orbit (in cyan). A 
trajectory on the chaotic attractor is depicted in Fig. 3(b), and 
the clearly irregular time series of the firing rate is shown in 
orange in Fig. 3(d). Using the same parameters we performed 
numerical simulations of the QIF neurons (l)-(2), finding a 
similar attractor and irregular dynamics as in the rate model, 
see Figs. 3(c) and (d). To obtain the time series shown in 
Figs. 3(d), we ran the QIF neurons numerically and, after a 
long transient, at time t = 0 , the rate model ( 12 ) was initiated 
with the values of r and v obtained from the population of 
QIF neurons. The time series of the two systems, which are 
initially close, rapidly diverge reflecting the chaotic nature of 


the system. 

Finally, to illustrate this chaotic state at the microscopic 
level. Fig. 3(e) shows the raster plot for 300 randomly chosen 
neurons, corresponding to the time series in Fig. 3(d). The 
irregular firing of neurons in Fig. 3(e) has some underlying 
structure that may be visualized ordering the same set of neu¬ 
rons according to their intrinsic currents as r/k < rjk+i, with 
k £ [1, 300], see Fig. 3(f). Clearly, the maxima of the firing 
rate coincide with the synchronous firings of clusters of neu¬ 
rons with similar 77 values. The size of these clusters is highly 
irregular in time, in concomitance with the chaotic behavior. 


E. Validity of the Lorentzian ansatz 


Thus far, we have shown that the LA (5) solves the conti¬ 
nuity equation (4), and confirmed that these solutions agree 
with the numerical simulations of the original system of QIF 
neurons (l)-(2). Here we further clarify why the LA holds for 
ensembles of QIF neurons. 

Transforming the voltage of the QIF neuron into a phase via 
Vj = tan(9j/2), the system (l)-( 2 ) becomes an ensemble of 
‘theta neurons’ [40]; 

Oj = (1 — cos Oj) + (1 + cos Oj)[rjj + Js(t) + /(f)]. (13) 


In the new phase variable 0 £ [0, 27r), the LA (5) becomes: 


P(8\v,t) = ^Re 


1 + a(? 7 ,f)e 


i6 " 


1 — a(j) 7 t)e 


w 


(14) 


where the function a(r/ 7 1) is related with 777 ( 77 , 1 ) as 




1 ~ w{rj,t) 

1 +w{rj,t)' 


(15) 


Equations (5) and (14) are two representations of the so- 
called Poisson kernel on the half-plane and on the unit disk, 
respectively. These representations are related via equation 
(15), that establishes a conformal mapping from the half-plane 
Re(w) > 0 onto the unit disk |a| < 1. In the next section we 
show that variables r and v can be related, via the same con¬ 
formal map, with the Kuramoto order parameter, which is a 
macroscopic measure of phase coherence [45, 46]. 

The key observation supporting the applicability of the LA 
is the fact that equation (14) turns out to be the ansatz dis¬ 
covered in 2008 by Ott and Antonsen [37]. According to the 
Ott-Antonsen (OA) theory, in the thermodynamic limit the dy¬ 
namics of the class of systems 

d t 0(r], t ) = ^(r?, f) + Irn [i/(?y, t)e~ z8 ] , (16) 


generally converges to the OA manifold (14). In our case, for 
equation (13), we have 0 ( 77 , f) = 1+rj+Js+I and H{r],t) = 
i (—1 + 77 + Js + /). Thus far, the convergence of (16) to 
the OA manifold has been proven only for H(i 7 , f) = H(t). 
This includes the well-known Kuramoto and Winfree models 
[47, 48], but not system (13). However, there are theoretical 
arguments [49] that strongly suggest that the OA manifold is 
also attracting for 77 -dependent H, as numerically confirmed 
by a number of recent papers using theta neurons [50-52] and 
other phase-oscillator models [53-57]. 
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FIG. 4. (color online). The conformal map (18) transforms the right 
half-plane onto the unit disk. This transformation and its inverse 
define a one-to-one mapping between the Kuramoto order parameter 
Z and the macroscopic quantities of the population of QIF neurons: 
r: firing rate, and v: mean membrane potential. Note that in the limit 
Z = e 1 ' 1 ' (full coherence, p(6) = S(9 — IF)) we recover v = i( 1 — 
e t9 )/(l+e l9 ) = tan('F/2), as the original mapping V =tan(0/2) 
between QIF- and theta-neurons dictates. 


F. Firing rate and Kuramoto order parameter 


There exists a mapping between the macroscopic variables 
r and v and the so-called Kuramoto order parameter Z. Equa¬ 
tion (15) relates, for each value of r), the firing rate and the 
mean membrane potential, both contained in w, with the uni¬ 
formity of the phase density, measured by a —note that a — 0 
in (14) yields a perfectly uniform density. The Kuramoto or¬ 
der parameter is obtained by integrating a over the whole pop¬ 
ulation as 

/ oo r2n r 

g(v) / p(0\r),t)e l9 d6dri 

-oo J 0 

/ oo 

g(ri)a*(r),t)dr) : (17) 

-OO 

where we assumed that the density p is the OA ansatz Eq. (14). 
The quantity Z can be seen as the center of mass of a popula¬ 
tion of phases distributed across the unit circle e . 

For a Lorentzian distribution of currents, g(rj), we may de¬ 
rive the exact formula (see Appendix B) 


Z = 


1 - W* 

1 + w* ’ 


( 18 ) 


relating the Kuramoto order parameter Z with the firing-rate 
quantity W = nr + iv. Figure 4 illustrates this conformal 
mapping of the right half-plane (r > 0) onto the unit disk 

(\Z\ < 1). 


m. CONCLUSIONS 

We have presented a method for deriving firing rate equa¬ 
tions for a network of heterogeneous QIF neurons, which is 
exact in the thermodynamic limit. To our knowledge, the re¬ 
sulting system of ordinary differential equations (12) repre¬ 
sents the first exact FREs of a network of spiking neurons. 


We would like to emphasize that the derivation of the FREs 
does not rely on assumptions of weak coupling, separation of 
time-scales, averaging, or any other approximation. Rather, 
the assumptions underlying the validity of equation (12) are 
that z) the QIF neurons be all-to-all connected, ii) heterogene¬ 
ity be quenched, and in) inputs be distributed according to a 
Lorentzian distribution. 

The last two assumptions may seem particularly restric¬ 
tive. Concerning (in), it must be stressed that for arbitrary 
distributions of quenched heterogeneity the LA (5) remains 
generally valid, and therefore equation (10) can be used to 
discern the stability of the network states (see Appendix C). 
Furthermore, in Appendix C we also show numerical simu¬ 
lations using uniform- and Gaussian-distributed inputs, that 
reveal very similar macroscopic dynamics in response to time 
varying inputs. Even relaxing assumption (ii), numerical sim¬ 
ulations with identical QIF neurons driven by external inde¬ 
pendent Gaussian noise sources show qualitative agreement 
with the FREs (12) (see Appendix D). In sum, the choice of 
a quenched Lorentzian distribution g(ij) is thus a mere math¬ 
ematical convenience, whereas the insights gained from the 
resulting firing-rate description are valid more generally. 

Our derivation represents a sharp departure from, and a ma¬ 
jor advance compared to previous studies in several regards. 
Firstly, in the past it has only been possible to calculate the 
approximate firing rate of networks of spiking neurons for 
stationary states, or for weak deviations of such states [9- 
12, 15, 20, 22]. The equations which result from these deriva¬ 
tions are difficult to solve, and typically require special nu¬ 
merical methods. In contrast, FREs (12) can be easily ana¬ 
lyzed and simulated, and exactly reproduce the behavior of the 
spiking network far from any fixed point and for arbitrary ex¬ 
ternal currents. Secondly, firing-rate descriptions traditionally 
assume that the activity of a population of neurons is equiva¬ 
lent to a set of uncorrelated stochastic processes with a given 
rate, see e.g. [13, 14]. However, in simulations of spiking net¬ 
works it is well known that the response of the network to 
non-stationary inputs generally involves some degree of spike 
synchrony. Recent theoretical work has sought to improve 
on classical rate models, which act as a low-pass filter of in¬ 
puts, by fitting them with linear filters extracted from the cor¬ 
responding Fokker-Planck equation for the network [20, 22]. 
These filters tend to generate damped oscillations when the 
external noise level is not too high, reflecting the presence of 
spike synchrony [15, 22], The FREs (12) also capture this 
phenomenon, and reveal that the underlying mechanism is an 
interplay between firing rate and subthreshold voltage. Fur¬ 
thermore, because these equations are exact, we can explore 
the full nonlinear response of the network, such as the gen¬ 
eration of chaotic states of synchronous bursting as shown in 
Fig. 3. 

Recently it has been shown that, not only Kuramoto-like 
models, but a much wider class of phase-models, evolve in 
the OA manifold defined by (14). Specifically, networks of 
pulse-coupled oscillators [48] and theta-neurons [50-52] al¬ 
low for an exact, low-dimensional description in terms of the 
Kuramoto order parameter (17). Although these later works 
use a finite-width phase coupling function that differs from the 
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synaptic coupling (3), the obtained low-dimensional descrip¬ 
tion in terms of the Kuramoto order parameter is analogous to 
ours, but in a different space [58], Indeed, we showed that the 
OA ansatz (14) is related, via the nonlinear transformation of 
variables (15), to the LA ansatz (5). Remarkably, this transfor¬ 
mation establishes an exact correspondence between the Ku¬ 
ramoto order parameter and a novel, biophysic ally maningful 
macroscopic observable which describes the firing rate and 
mean membrane potential of the neuronal network. Interest¬ 
ingly, the low-dimensional description in terms of firing rates 
seems to be a more natural description for networks of spik¬ 
ing neurons, compared to that in terms of the Kuramoto order 
parameter. The firing rate equations (12) take a surprisingly 
simple form in the LA manifold, what makes them a valuable 
tool to explore and and understand the mechanisms governing 
the macroscopic dynamics of neuronal networks. 

Finally, since the OA ansatz is the asymptotic solution for 
systems of the form in equation (16), applying the change of 
variables V = tan(0/2) automatically implies that the LA 
should hold for populations governed by: 

Vj = A + BinjrfVj + C(r,j,t) V?, (19) 


where A, B, and C are related with 12 and H as A = 
[12 + Im(iT)]/2, B = —Re(iT), C = [12 - Im(iT)]/2. No¬ 
tably, equation (19) defines a wide family of ensembles of QIF 
neurons. 


Therefore, the LA is actually valid for more general net¬ 
works beyond the one investigated here — rjj in Eq. (19) may 
also be a vector containing several forms of disorder. As par¬ 
ticularly relevant cases, in Appendix E we provide the FREs 
governing the dynamics of an excitatory network with both, 
heterogeneous inputs ?) and synaptic weights J, as well as pair 
of interacting excitatory and inhibitory populations of QIF 
neurons. In addition, according to Eq. (19), the LA is also 
valid if synapses are modeled as conductances, in which case 
the reversal potentials may be distributed as well. Moreover, 
the role of gap-junctions or synaptic kinetics may be consid¬ 
ered within the same framework. 
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APPENDIX A: NUMERICAL SIMULATIONS 

To numerically simulate the population of QIF neurons (1) 
we used the Euler method, with time step dt = 10 -4 . The 
population had N = 10 4 neurons, and the Lorentzian dis¬ 
tribution (11) was deterministically generated using: ijj = 
?7+Aarctan[7r/2(2j —iV—1)/(7V+1)], where j = 
and A = 1. 

The time it takes for the membrane potential of a QIF neu¬ 
ron (with I 3 > 0) to reach infinity from a given positive 
value of the membrane potential is ai'ctan( /V p )/^/Tj. 
For yJTj <C V p , this expression can be approximated as: 

arctan ( /V p ) 1 

7 % “ V 

In simulations we considered V p = —V r = 100, and used 
the previous approximation. Thus, the time for the neurons to 
reach infinity from V p is 1/V P w 10~ 2 . Additionally, the time 
taken from minus infinity to V r is 1/V r 10 -2 . Numerically, 
once a neuron’s membrane potential Vj satisfies Vj > V p the 
neuron is reset to — Vj and held there for a refractory time 
2 /Vj. The neurons produce a spike when Vj reaches infinity, 
i.e. a time interval 1 /Vj after crossing V p . The exact time of 
the spike can not be evaluated exactly in numerical simula¬ 
tions, since dt is finite. However, simulations agree very well 
with the theory provided that 1 /V p dt. 

To evaluate the mean membrane potential v, the popula¬ 
tion average was computed discarding those neurons in the 
refractory period. The choice V p = —V r is not needed, but 
in this way the observed mean membrane potential v agrees 
with the theory, consistent with the definition in (8). The mean 
synaptic activation (3) was evaluated using the Heaviside step 
function a T (t) = 0(r — f)/r with r = 10~ 3 . The instan¬ 
taneous firing rates were obtained binning time and counting 
spikes within a sliding time window of size St = 2 x 1CH 2 
(Figs. 2, C2, C3, Dl) and St = Ax 10" 2 (Fig. 3). 


APPENDIX B: PROOF OF EQUATION (18) 


The inverse of Eq. (18) reads: 


W = 


1 - Z* 

1 + z* 


(Bl) 


Recalling Eqs. (7) and (9), we write the macroscopic field W 
as 


W(t) = 



w(v,t)g{v)dr)- 


where w = x + iy. Inserting the conformal mapping w = 
(1 — a)/(l + a) — the inverse of Eq. (15)- we get 


W{t) = 


/*oo 

'1 - a(r), t)~ 

I —00 

1 + a(r],t)_ 


g{v)dv■ 
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Using the geometric series formula, and grouping the powers 
of a, we may evaluate the integrals obtaining 

W{t) = 1- 2Z*(t) + 2Z%(t) - 2 z;(t) + • • • , (B2) 

where the Z m are the generalized order parameters [63]: 

/ oo 

g{il) / p(6\y,t)e lme d6 dg 

-oo J 0 

/ oo 

-OO 

Since, for Lorentzian 77 ( 77 ), Z m (t) = [a*(g = fj — = 

Z(t) m , we can revert the power series in (B2) and obtain the 
result in (Bl). 


APPENDIX C: RESULTS FOR ARBITRARY 
DISTRIBUTIONS OF CURRENTS 


In this Appendix we compare the results in the main text, 
obtained using a Lorentzian distribution of currents < 7 ( 77 ), with 
other distributions. The results are qualitatively similar as ev¬ 
idenced by Figs. Cl, C2 and C3 (cf. Figs. 1, 2, and 3 in the 
main text). 

Note that if a particular distribution 77 ( 77 ) has 2 n poles (all of 
them off the real axis) one can readily obtain the FREs con¬ 
sisting of n complex-valued ordinary differential equations. 
Even if 7 /( 77 ) does not fulfill this condition, as in the case of 
Gaussian or uniform distributions, the Lorentzian Ansatz still 
holds (see below). 

Moreover, it is possible to efficiently simulate the dynamics 
of a population by integrating Eq. (10) for a sample of 77 val¬ 
ues that approximate a particular 77 ( 77 ) [64]. However, if 77 ( 77 ) 
is a discrete distribution (or has a discrete part), —as in the 
case of the Dirac delta function 77 ( 77 ) = (>(77 — fj) —, the so- 
called Watanabe-Strogatz theory applies [49, 65, 66 ], and the 
‘Lorentzian manifold’ is no longer attracting since the dynam¬ 
ics is extremely degenerate, akin to a Hamiltonian system. 


Steady states 

Considering I(t ) =0 in Eqs. (1) and (2), it is possible to 
find the solutions with steady firing rate r(t) = ro for arbitrary 
distributions of currents. 

As we pointed out in the main text, the stationary solution 
of the continuity equation (4) for those neurons that are intrin¬ 
sically active (77 > —Jtq) must be is inversely proportional to 
their speed, that is 


Po(V\tj) 


c(il) 

V 2 + 77 + Jr 0 ’ 


where the normalizing constant is 77 ( 77 ) =■ y/g + Jr^/it. The 
firing rate for each value of 77 is then 

ro(v) = Po(V —f oo|? 7 )U(V -A 00 ( 77 , f) = c(v), 


and integration over all the firing neurons, gives the self- 
consistent condition for r 0 : 


ro 



c(v)g(i 7 ) di 7 


(Cl) 


which is valid for any distribution 77 ( 77 ) (an analogous expres¬ 
sion was obtained in Eq. (45) of [42]). For Lorentzian 77 ( 77 ), 
the integral in (Cl) can be evaluated, and solving the result¬ 
ing equation for ro one obtains the steady states in agreement 
with the result of setting f = v = 0 in the FREs (12). 


Linear stability analysis of steady states 


Assuming that the LA captures the actual dynamics of the 
population of QIF neurons, one can also investigate the stabil¬ 
ity of steady states for arbitrary distributions 77 ( 77 ). Indeed, the 
evolution of infinitesimal perturbations from a steady state 


7770(77) = £0(77) +iyo{g) 


\J 77 + Jr 0 if 77 > — Jr 0 
—iy/—g — Jr 0 if 77 < —Jr 0 


is determined by the linearization of Eq. (10): 


dt 5 w(y, t) = i(JSr — 2wq6w) (C2) 

where Sr(t) = ( 27 T )” 1 [75x77(77, i) + Sw*(g,t)]g(g)dg. 

To find the discrete spectrum of eigenvalues [67], we use 
the ansatz 5 w(g, t) = b(jf)e xt in (C2) and obtain 


[A+2?T;o(77)i]6(77)e At 


• T POO 

l — J [b(r])e xt +b*{y)e x * t }g('n)d 7 i. 


Solving the equation for be xt , we find 


b( g)e xt 


1 J 
27770(77) — iA 27 t 



[b(r])e xt +b* (g)e x * ^g (g)dy. 


After summing the complex conjugate at each side of the 
equation, we multiply by 77 ( 77 ) and integrate over 77 . This 
allows one to cancel out the integrals J(be xt + b*e x *)gdg. 
Then, imposing self-consistency, we find 



1 

2tt7o(?7) — *A 


+ c.c. 


g(v)dg. 


(C3) 


Local bifurcations are associated with marginal stability of the 
fixed points: Re(A) —7 0. In correspondence with the results 
obtained for a Lorentzian distribution in the main text, we ex¬ 
pect bifurcations of saddle-node type associated with A = 0. 
Therefore, from Eq. (C3) we find that the loci of the saddle- 
node bifurcations are located at 


JSN 


2 n 

"OO x 0 (y) 

—00 (uioG)F 


g(v)dg 


(C4) 


Since xo = V 7 ? + Jr 0 for 77 > — Jr 0 , and zero otherwise, 
we can restrict the integration over 77 to the range (—J 7 ’o, 00 ). 
Now multiplying Eq. (C4) by Eq. (Cl) we get, after defining 
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FIG. Cl. (color online). Phase diagrams for (a) uniform and (b) 
Gaussian distributions of currents g(ji). The cyan-shaded region 
represents the bistable region, with the solid lines corresponding to 
saddle-node bifurcations analytically obtained from Eqs. (C5) and 
(C 6 ). Square symbols are estimations of the bifurcations loci ob¬ 
tained from by direct numerical simulations of N = 10 4 QIF neu¬ 
rons. Triangle symbols indicate parameter values used in numerical 
simulations of Figs. C2 and C3. 


£ = JsnI" o- t wo equations that permit to find the locus of the 
saddle-node bifurcations systematically for arbitrary distribu¬ 
tions of currents: 


JSN = 


2 tt 


J^°ri- 1/2 g(V-0dv 


with £ obtained solving 


£ = 


2 / 0 oo? ? 1/2 g( 7 ?-£)^ 

fo°° v~ 1/2 g(v - £)<*?' 


(C5) 


(C 6 ) 


Bistability region for uniform and Gaussian distributions 


Equations (C5) and (C 6 ) are particularly easy to solve for 
uniform distributions of currents: 


g(v) = 



for I 77 — r)| < 7 
otherwise 


After defining the rescaled parameters fj = 77/7 and J = 
J/y/ 7, we find two branches of saddle-node bifurcations em¬ 
anating from the cusp point at fj = —1/3: 

f( 1 ) _ 2n 
SN ~ V / 3^T3 

^_ 

yJfj+l + 2^±+fj 2 - \jf)-l + 2 s J\+rj 2 

These functions are plotted in Fig. SI (a). Numerical simula¬ 
tions using the network of QIF neurons confirm the correct¬ 
ness of these boundaries (see the square symbols in Fig. Cl). 
For Gaussian distributions, 


g 0?) = 


1 r -(v-v) 2 /(^ 2 ) 

V^a 


the solutions of (C5) and (C 6 ) can be numerically found. The 
results are shown in Fig. Cl(b). Again, there is a perfect 
agreement between Eqs. (C5) and (C 6 ) —derived assuming 
the Lorentzian Ansatz— and the numerical estimations ob¬ 
tained simulating the network of QIF neurons. 


Excitatory networks with external periodic currents 

Firing rate equations (12) predict the existence of a stable 
focus in the shaded region of Fig. 1. Trajectories attracted to 
this fixed point display oscillations in the firing rate and mean 
membrane potential due to the transient synchronous firing of 
the QIF neurons. 

When an external periodic current is injected to all neurons 
in the network, the spiral dynamics around the fixed point is 
responsible for the bursting behavior observed in Fig. 2, as 
well as for the emergence of the macroscopic chaos shown in 
Fig. 3. Here we investigate whether similar phenomena oc¬ 
cur when an external periodic current of the same frequency 
is injected in an excitatory network with either uniform or 
Gaussian-distributed currents. 

We introduced a sinusoidal forcing I (t) = Iq sin {tot) ob¬ 
serving a behavior qualitatively identical to the one reported 
in the main text with Lorentzian < 7 ( 77 ). Under a low-frequency 
forcing the system displays periodic bursting, see Fig. C2, 
provided the parameters are set inside the bistable region (see 
the black triangles in Fig. Cl). Furthermore, note that the 
range of firing rates and mean membrane potentials in Fig¬ 
ure C2 are similar to those of Fig. 2 —for Lorentzian distribu¬ 
tions of currents. 

As for the simulations shown in Figure (3), we next increase 
the frequency of the injected current up to w = 7 r to investi¬ 
gate the emergence of macroscopic chaos in a network with 
Gaussian-distributed currents. Fig. C3 shows the emergence 
of an apparently chaotic state. The observed dynamics is sim¬ 
ilar to that of Fig. 3, which suggests it is chaotic. This type 
of chaos persists in the thermodynamic limit N —> 00. On 
top of this, highly-dimensional but weakly chaotic dynamics 
is probably also present due to ‘residual’ finite-size effects. 


APPENDIX D: IDENTICAL QIF NEURONS DRIVEN BY 
INDEPENDENT NOISE TERMS 


We compare now the results obtained above and in the main 
text for quenched heterogeneity < 7 ( 77 ) with the results for iden¬ 
tical neurons < 7 ( 77 ) = 6(17 — fj) driven by noise. Specifically, 
the inputs currents are now taken as 

Ij = V+ Jr(t) +£j(f), (Dl) 

where £j are independent white noise terms with expected val¬ 
ues (£j(f)) = 0 , and (£y(f)£fc(f')) = 2D s jk S{t - t'). 

In the thermodynamic limit, the density p(V, t) obeys the 
Fokker-Planck equation: 

dtp + dy [{V 2 +fj + Jr)p\ = D dyp. (D2) 

The Lorentzian ansatz is not a solution of this equation, 
but we demonstrate here that the phenomena observed with 
quenched heterogeneity arise also with independent noise 
sources. This qualitative similarity at the macroscopic level 
between quenched Lorentzian heterogeneity and Gaussian 
noise has been noted in previous work [60]. 
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FIG. C2. Numerical simulations of the excitatory network of QIF neurons, Eqs. (1), and (2). with (a-d) uniform and (e-h) Gaussian distributions 
of currents. As in Fig. 2, at time t = 0, an external sinusoidal current I ( t ) = Iq sin(uT) —shown in panels (d,h)— is applied to all neurons. 
Panels (a,e) show the time series of the firing rate, and (b.f) the mean membrane potential of all cells. In panels (c,g) we depicted the raster 
plots of 300 randomly chosen neurons. Parameters correspond to the black triangle symbols in Fig. Cl: (a-d) fj = — 1, J = 8, 7 = 1, Jo = 1; 
(e-h) 77 = -2, J = 10, o = 1, J 0 = 1.5. 




FIG. C3. Chaotic state in a network of N = 10 4 QIF neurons with 
Gaussian-distributed currents. Network parameters correspond to the 
yellow-triangle symbol in the phase diagram of Fig. S1 (b). Like in 
Fig. 3, neurons receive a common periodic current I(t ) = Jo sin(cuf) 
of frequency w = 7r. Parameters: fj = —1, J = 8, a = 1, Jo = 1.5. 


to each fixed point (see figure caption for details). 

Additionally, in order to investigate and compare the dy¬ 
namical behavior of the identical noise-driven neurons with 
that of the FREs (12), an external current of intensity Iq = 3 
is injected to all neurons at time t = 0, like in Fig. 2. In 
Fig. Dl(b,c) the time series of the firing rate and the mean 
membrane potential clearly display damped oscillations after 
the injection of the current, confirming the existence of a sta¬ 
ble focus, exactly as observed in the FREs (12). The existence 
of a stable focus reflects the presence of transient spike syn¬ 
chrony in the network, as seen by the raster plot in Fig. Dl(d). 
Remarkably, the raster plot and the time series closely resem¬ 
ble those of Fig. 2. The resemblance is not just qualitative, 
but rather there is a near quantitative fit between the network 
of QIF neurons driven by Gaussian noise and the FREs (12), 
which were derived assuming quenched Lorentzian hetero¬ 
geneity (orange curves). 


APPENDIX E: MODEL GENERALIZATIONS 


To illustrate the potential of the LA for investigating more 
sophisticated networks, here we provide the low-dimensional 
FREs corresponding to an excitatory network of QIF neurons 
with independently distributed currents and synaptic weights, 
and to two interacting populations of excitatory and inhibitory 
QIF neurons with distributed currents. 


A subtle point in Eq. (D2) is that its nondimensionalization 
entails a different (compared to A) rescaling with D: V = 
V/D 4 / 3 , fj = rj/D 2/3 , J = J/D 4 / 3 , t = tD 1/3 (implying 
r = r/D - 1 / 3 ). 

Numerical simulations reveal the existence of a region of 
bistability, see Fig. Dl(a), analogous to the one observed for 
quenched heterogeneity, cf. Figs. 1(a) and Cl. Obviously, true 
bistability only holds in the thermodynamic limit, while what 
we observe are rather exceedingly long residence times close 


Excitatory population with heterogeneous currents and synaptic 

weights 

As a first example, let us assume that both the currents 77 
and the synaptic weights J are distributed —this type of het¬ 
erogeneity was also considered in [42]. The input currents 
read then [ 68 ] 

=Vj +Jjr(t) +/(t). 
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(b) 


(c) 


<d) 


(e) 



time 


FIG. Dl. (color online). Numerically obtained phase diagram (a) and time series (b-d) for a network of identical QIF neurons driven by 
independent white noise terms. The square symbols in panel (a) represent the boundaries that enclose the region of bistability between a 
high and a low activity state. The boundaries have been obtained using a network of 10 4 neurons, and by numerical continuation of a low 
activity solution (filled squares) and a high activity solution (empty squares). Specifically, using a noise intensity D = 1 and a particular 
value of J, the system was initialized either at fj = —5 or at fj = 0 and, after t = 10 time units, parameter fj was decreased/increased an 
amount 0.025/0.1, respectively. This continuation was made until the relative change of two successive values of the averaged firing rate 
(time-averaged over the last time unit for each parameter value) was larger than 50%. Other parameters for the numerical simulations are 
the same as in all other figures, and are described in the Material and Methods section. The triangle symbol indicates the parameter value, 
{fj/ D 2 ^ 3 , J/D 4 / 3 ) = (—5,15), corresponding to the numerical simulations shown in black in the right panels —for the simulations we used 
D = 1. To facilitate the comparison with the FREs (12), the orange curves show the time series of the FREs, using the same parameters 
(■ fj , J) = (—5,15), but with A = 1 


For simplicity, let us additionally assume that rj and J are 
distributed independently, with a joint distribution p(f), J) = 
g{rj)h{J). In the simplest situation of Lorentzian g{rj) and 
h(J) 

t t _ A / ?r Uf t\ _ V /' K 

(v-fj ) 2 + a 2 ’ h ^ J) (J-J ) 2 + r 2 ' 


the problem is extremely simplified using the Lorentzian 
ansatz, which now trivially reads: 




tt[V - 2 /( 77 , J, t )] 2 + x{rj, J, t ) 2 ' 


The firing rate and mean membrane potential are determined 
only by the value of w = x + iy at the poles in the lower half 
planes: 



FIG. El. (color online). Phase diagram for an excitatory population 
with heterogeneous currents and synaptic weights obtained using the 
FREs (El). 


r{t) + iv(t ) 


= JJ w{r],J,t)g{g)h{J)drjdJ 


= iv{rj — iA, J — iT, t). 


Firing-rate equations for a pair of Excitatory-Inhibitory 
populations 


Finally, evaluating Eq. (10) in the main text at the poles (77 = 
77 — iA, J = J — tT) we get the exact FREs: 

r = A/- 7 T + Tr/ir + 2rv, (Ela) 

v = v 2 + fj + Jr + I(t) — 7 r 2 r 2 . (Elb) 

These equations contain simply an extra term ‘+rr/ 7 r’ com¬ 
pared to equations (12). Figure El shows the bistability 
boundaries obtained from (El) for different ratios of 1 and 
A, and I{t) = 0. Note that the region of bistability shifts 
to lower values of fj/A and to higher values of J/v/A as the 
level of heterogeneity in the synaptic coupling T is increased. 


The microscopic state of each population of QIF neurons is 
characterized by the membrane potentials {TA ( ' e ’ z - ) }j=i 
which obey the following ordinary differential equations: 



(y . ( ed ) ) 2 + 7 je,i ) 


if y\ (e,i) > V/.thcn vj e,i) <- V r . 


(e i) 

Here, Vj ’ represents the membrane potential of neuron j in 
either the excitatory (e), or in the inhibitory population (i). 
The external currents for the excitatory and inhibitory popula- 
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tions are, respectively: 

/j 6) = Vf + JeeS {B \t) ~ JieS®(t) + 

if = Tjf + Jei/ e \t) - + / (i) (i), 

where the synaptic weights are J ee , Jj, ; , J ie , J ei . Finally, the 
mean synaptic activation for each population is 

j =1 

Here, (^)( e >») is the time of the fcth spike of yth neuron in 
either the excitatory (e), or in the inhibitory population (*). 
Additionally, 5{t) is the Dirac delta function, and a T (t) is the 
normalized synaptic activation caused by a single pre-synaptic 
spike with time scale r, e.g. a T (t) = e _t / T /r. 


dt'a r (t-t')S(t'-(^) (eA ). 


It is straightforward to apply the LA and the method de¬ 
scribed in the main text, to obtain the firing rate equations 
corresponding to the two-population model. Considering the 
limit of infinitely fast synapses, r —> 0, we get s( e >®) (£) = 
r (e,0(£) Finally, assuming the Lorentzian distributions of 
currents for both populations: 


9e,i(v) 


1 A e>< 

7T (t? - rje,i ) 2 + A 2 e i ’ 


we obtain the firing-rate equations: 

f( e ) = Ag/ 7 r + 2 r^v^ e \ 

v {e) = (f (e) ) 2 + % + Jeer^ e) - J le r« + 
fW = Aj/7r + 

= (r ,«) 2 + fji + J ei rW + J«(t) - (nr^) 2 . 



